﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace TopRelax
{
    public class SORelax
    {
        //код для простой итерации
        public double[,] u; //искомое решение
        public double[,] b;
        public int n; //число разбиений по икс
        public int m; //число разбиений по игрек
        public double h; // шаг по икс
        public double k; // шаг по игрек
        
        public double w; // параметр метода

        private double di;

        public SORelax()
        {
            n = m = 0;
            h = k = 0;
            setDi();
            u = null;
            b = null;
            
            w = 0;
        }
        public SORelax(int n, double h, int m, double k, double[,] b, double[,] u)
        {
            this.n = n;
            this.m = m;
            this.h = h;
            this.k = k;
            setDi();
            this.u = new double[n + 1, m + 1];
            this.b = new double[(n + 1) , (m + 1)];
            for (int i = 0; i <= n; i++)
            {
                for (int j = 0; j <= m; j++)
                {
                    this.u[i, j] = u[i, j];
                    this.b[i,j] = b[i,j];
                }
            }
            
            setW();

        }
        private void setDi()
        {
            di = -2.0 * (1.0 / (h * h) + 1.0 / (k * k));
        }
        private void setW()
        {
            double[,] matr = new double[(n-1) * (m-1), (n-1) * (m-1)];
            double d1 = 1.0 / (h * h);
            double dn = 1.0 / (k * k);

            //заполнение матрицы
            for (int c = 0; c < (n-1)*(m-1) ; c++)
            {
                matr[c, c] = 0;
                if (c != 0 && ((c+1)%(n-1) - 1 !=0)) matr[c, c - 1] = d1/di;
                if (c != (n - 1) * (m - 1) - 1 && ((c + 1) % (n - 1) - 1 != -1)) matr[c, c + 1] = d1/di;
                if (c > n - 2) matr[c, c - (n - 1)] = dn/di;
                if (c <= (n - 1) * (m - 1) - n) matr[c, c + (n - 1)] = dn/di;
            }

            //поиск собственных чисел
            double[] wr =new double[(n-1)*(m-1)];
            double[] wi;
            double[,] wl,wu;
            alglib.rmatrixevd(matr, (n - 1) * (m - 1), 0, out wr, out wi,out wl,out wu);

            //поиск максимума
            int Imax = 0;
            for (int i = 0; i < (n - 1) * (m - 1); i++)
            {
                if (Math.Abs(wr[Imax]) < Math.Abs(wr[i])) Imax = i;   
            }

            //формула оптимального параметра
            w = 2.0 / (1 + Math.Sqrt(1 - wr[Imax] * wr[Imax]));
        }

        private double findMaxAbs(double[,] d1, double[,] d2, out int Imax, out int Jmax)
        {
            Imax = 1;
            Jmax = 1;
            for (int i = 1; i < n; i++)
                for (int j = 1; j < m; j++)
                    if(Math.Abs(d1[i,j] - d2[i,j])>Math.Abs(d1[Imax,Jmax] - d2[Imax,Jmax]))
                    {
                        Imax = i;
                        Jmax = j;
                    }
                
            return Math.Abs(d1[Imax,Jmax] - d2[Imax,Jmax]);
        }

        public int method(double e,out double ereal, int count,int Imax,int Jmax)
        {
            
            int s = 0;
            double[,] oldu = new double[n+1,m+1];
            
            double dl = 10;//точности метода
            double d1 = 1.0 / (h * h);
            double dn = 1.0 / (k * k);
            
            while (dl > e && s < count)
            {
                s++;                  
                
                for (int p = 0; p <= n; p++)
                    for (int j = 0; j <= m; j++)
                        oldu[p, j] = u[p, j];

                for (int j = 1; j < m; j++)
                {
                    for (int p = 1; p < n; p++)
                    {
                        u[p, j] = ( w*( -d1*(u[p-1,j]+oldu[p+1,j]) -dn*(u[p,j-1]+oldu[p,j+1]) + b[p,j]) )/di + (1-w)*oldu[p,j];
                    }
                }
            
                dl = findMaxAbs(u, oldu, out Imax, out Jmax);
            }
            ereal = dl;
            return s;
        }


    }
}
